%% BioMed_Central_Tex_Template_v1.06
%%                                      %
%  bmc_article.tex            ver: 1.06 %
%                                       %

%%IMPORTANT: do not delete the first line of this template
%%It must be present to enable the BMC Submission system to
%%recognise this template!!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                     %%
%%  LaTeX template for BioMed Central  %%
%%     journal article submissions     %%
%%                                     %%
%%         <14 August 2007>            %%
%%                                     %%
%%                                     %%
%% Uses:                               %%
%% cite.sty, url.sty, bmc_article.cls  %%
%% ifthen.sty. multicol.sty		       %%
%%				      	               %%
%%                                     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                 %%	
%% For instructions on how to fill out this Tex template           %%
%% document please refer to Readme.pdf and the instructions for    %%
%% authors page on the biomed central website                      %%
%% http://www.biomedcentral.com/info/authors/                      %%
%%                                                                 %%
%% Please do not use \input{...} to include other tex files.       %%
%% Submit your LaTeX manuscript as one .tex document.              %%
%%                                                                 %%
%% All additional figures and files should be attached             %%
%% separately and not embedded in the \TeX\ document itself.       %%
%%                                                                 %%
%% BioMed Central currently use the MikTex distribution of         %%
%% TeX for Windows) of TeX and LaTeX.  This is available from      %%
%% http://www.miktex.org                                           %%
%%                                                                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\NeedsTeXFormat{LaTeX2e}[1995/12/01]
\documentclass[10pt]{bmc_article}



% Load packages
\usepackage{cite} % Make references as [1-4], not [1,2,3,4]
\usepackage{url}  % Formatting web addresses
\usepackage{ifthen}  % Conditional
\usepackage{multicol}   %Columns
\usepackage[utf8]{inputenc} %unicode support
%\usepackage[applemac]{inputenc} %applemac support if unicode package fails
%\usepackage[latin1]{inputenc} %UNIX support if unicode package fails
\urlstyle{rm}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
%%                                             %%
%%  If you wish to display your graphics for   %%
%%  your own use using includegraphic or       %%
%%  includegraphics, then comment out the      %%
%%  following two lines of code.               %%
%%  NB: These line *must* be included when     %%
%%  submitting to BMC.                         %%
%%  All figure files must be submitted as      %%
%%  separate graphics through the BMC          %%
%%  submission process, not included in the    %%
%%  submitted article.                         %%
%%                                             %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\def\includegraphic{}
\def\includegraphics{}



\setlength{\topmargin}{0.0cm}
\setlength{\textheight}{21.5cm}
\setlength{\oddsidemargin}{0cm}
\setlength{\textwidth}{16.5cm}
\setlength{\columnsep}{0.6cm}

\newboolean{publ}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                              %%
%% You may change the following style settings  %%
%% Should you wish to format your article       %%
%% in a publication style for printing out and  %%
%% sharing with colleagues, but ensure that     %%
%% before submitting to BMC that the style is   %%
%% returned to the Review style setting.        %%
%%                                              %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%Review style settings
\newenvironment{bmcformat}{\begin{raggedright}\baselineskip20pt\sloppy\setboolean{publ}{false}}{\end{raggedright}\baselineskip20pt\sloppy}

%Publication style settings
%\newenvironment{bmcformat}{\fussy\setboolean{publ}{true}}{\fussy}



% Begin ...
\begin{document}
\begin{bmcformat}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the title of your article here     %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{LSPR: An integrated periodicity detection algorithm for unevenly sampled temporal microarray data}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors here                   %%
%%                                          %%
%% Ensure \and is entered between all but   %%
%% the last two authors. This will be       %%
%% replaced by a comma in the final article %%
%%                                          %%
%% Ensure there are no trailing spaces at   %%
%% the ends of the lines                    %%     	
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\author{Rendong Yang\equalcontributors$^1$%
       \email{Rendong Yang - cauyrd@gmail.com}%
      \and
         Chen Zhang\equalcontributors$^2$%
         \email{Chen Zhang - zcreation@yahoo.cn}
       and
         Zhen Su\correspondingauthor$^1$%
         \email{Zhen Su\correspondingauthor - zhensu@cau.edu.cn}%
      }


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors' addresses here        %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\address{%
    \iid(1)Division of Bioinformatics, State Key Laboratory of Plant Physiology and Biochemistry, College of Biological Sciences, China Agricultural University, Beijing 100193, China\\
    \iid(2)Department of Applied Mathematics, College of Science, China Agricultural University, Beijing 100083, China
}%

\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Abstract begins here                 %%
%%                                          %%
%% The Section headings here are those for  %%
%% a Research article submitted to a        %%
%% BMC-Series journal.                      %%
%%                                          %%
%% If your article is not of this type,     %%
%% then refer to the Instructions for       %%
%% authors on http://www.biomedcentral.com  %%
%% and change the section headings          %%
%% accordingly.                             %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{abstract}
        % Do not use inserted blank lines (ie \\) until main body of text.
        \paragraph*{Background:} Periodic processes are widespread in biology, such as circadian rhythm and cell cycle. Identifying genes with periodic expression patterns is the first step to study these processes. Most established periodicity inference algorithms can only be applied to evenly spaced time-course microarray data. However, as for limitations of technique and exorbitant budget on microarray experiments, the time points are relatively insufficient and unevenly spaced when missing point occurs. Therefore, the development of a new high-resolution algorithm is essential for analyzing irregularly sampled temporal expression profiles.

        \paragraph*{Results:} We proposed a three-step combining periodicity detection algorithm named LSPR. After data preprocessing, our method integrates Lomb-Scargle periodogram to estimate the periods for each transcript, and then employs harmonic regression to evaluate the significance of cyclic components in the time-series. Results from simulated unevenly sampled temporal data sets shows that LSPR is best-performed compared with two widely used detecting methods: COSOPT and Lomb-Scargle periodogram. Furthermore, we have applied our algorithm to \emph{Arabidopsis} diurnal expression data and Yeast cell cycle expression data and compare its performance with existing well-established algorithms. The results indicate that LSPR are capable of identifying periodic transcripts more precisely than previous reports, especially for short and unevenly spaced time-course expression data.

        \paragraph*{Conclusions:} Testing on well-defined simulating and real molecular data, our proposed method is efficient to identify sinusoidal and non-sinusoidal periodic patterns in short, noisy and unevenly sampled time-series. LSPR algorithm is implemented by MATLAB software and is available for download at \url{http://bioinformatics.cau.edu.cn/LSPR}
\end{abstract}



\ifthenelse{\boolean{publ}}{\begin{multicols}{2}}{}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Main Body begins here                %%
%%                                          %%
%% The Section headings here are those for  %%
%% a Research article submitted to a        %%
%% BMC-Series journal.                      %%
%%                                          %%
%% If your article is not of this type,     %%
%% then refer to the instructions for       %%
%% authors on:                              %%
%% http://www.biomedcentral.com/info/authors%%
%% and change the section headings          %%
%% accordingly.                             %%
%%                                          %%
%% See the Results and Discussion section   %%
%% for details on how to create sub-sections%%
%%                                          %%
%% use \cite{...} to cite references        %%
%%  \cite{koon} and                         %%
%%  \cite{oreg,khar,zvai,xjon,schn,pond}    %%
%%  \nocite{smith,marg,hunn,advi,koha,mouse}%%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%
%% Background %%
%%
\section*{Background}
 Time-course microarray data are usually generated to study periodic processes such as cell cycle and circadian rhythm. The design of temporal microarray experiment is usually with low sampling rate and short time-series length. A typical circadian expression profile has less than twenty time-points \cite{duffield2003dna}. The reasons why most common type of temporal data available is short time-series data are high costs of the arrays and dampening of the circadian rhythm \cite{bar2004analyzing}. After generating these time-series microarray data, specially designed statistical methods are proposed to identify genes with periodic patterns from the whole genome transcripts.

 Limited sampling accentuates the difficulties associated with classical time-series analyses \cite{wang2008short}. Challenges for analyzing such short time series are high dimensionality accompanied by a small sample size (on the order of $10^{4}$ gene IDs per microarray) \cite{braga2007fads} and the unavoidable noise has more influence on the analysis of short time-series than on long time-series data \cite{ernst2005clustering}.

 To overcome these challenges, a variety of algorithms have been developed to identify the transcripts with periodic patterns from temporal expression data \cite{Yamada2007}. Generally, most of these methods are proposed based on traditional spectral estimation technology, such as autocorrelation \cite{wijnen2005molecular}, Fourier transformation \cite{de2005comparison}, Fisher's g-test \cite{wichert2004identifying} and maximum entropy algorithm \cite{langmead2002maximum}. These methods have been adapted to analyze circadian expression data of Mouse, \emph{Drosophila} and cell cycle data of Yeast. Recently, an algorithm called ARSER was proposed based on harmonic regression and autoregressive spectral estimation \cite{yang2010analyzing}. By testing on simulation and \emph{Arabidopsis} microarray data, it is shown ARSER can successfully identify rhythmic transcripts from microarray data.

 Although all these spectral estimation methods have their own advantages, they generally required the input signals were equally sampled. To look for a cyclic expression profile with one of these methods, the data must be collected at equal intervals (e.g., 10 min for cell cycle data, 4 hrs for circadian rhythm data) for the whole cycles (e.g., 1 day, 2 days, etc.) \cite{de2006chronobiometry}. However, for some biological experiments, there can be extra values or missing time points in the time series.

 For irregularly sampled time series, ``interpolated'' values can be inserted into the missing data points to satisfy the requirement of spectral estimation methods \cite{wichert2004identifying}. However, this way may not be a reflection of the true values that have been obtained. Recently, specially designed spectral estimation methods for analyzing unevenly sampled time-series have been developed. Alan \emph{et al}. \cite{liew2007spectral} proposed a spectral estimation method based on signal reconstruction in a shift invariant signal space to analyze Plasmodium falciparum Yeast expression data. Ahdesm$\ddot{a}$ki \emph{et al}. \cite{ahdesm?ki2007robust} proposed the robust regression to detect periodicity in non-uniformly sampled time-course gene expression data. In addition, Lomb-Scargle periodogram has been applied to unevenly sampled gene expression data for spectral estimation \cite{glynn2006detecting}. It is a direct method to treat the missing values and irregularly sampled time series and has been used by Dequeant \emph{et al}. \cite{dequeant2006complex} to detect rhythmic transcripts of mouse segmentation clock expression profiles. A comparison study of periodicity detection methods for irregularly sampled data concluded that Lomb-Scargle method performed better than most existing methods \cite{wentao2008detecting}.


 Alternatively, some studies use least square technique to fit the time series by sinusoidal models. A widely used algorithm called COSOPT \cite{straume2004dna} is designed based on such principle and this method has been successfully applied to analyze circadian microarray data of \emph{Arabidopsis} \cite{Edwards2006}, \emph{Drosophila} \cite{ceriani2002genome}, and mammalian systems \cite{10.1371/journal.pcbi.0020016, Hughes2009}. The COSOPT method fits the time series by cosine waves with varying phases and periods, and measures the significance of goodness-of-fit. Previous algorithms to analyze unevenly sampled data are subject to noise and not powerful for short time series. There is a need for an efficient method to identify periodicity in unevenly sampled short time series.

 In this study, we proposed a new periodicity identification algorithm based on Lomb-Scargle periodogram and harmonic regression method for unevenly sampled time series. Our algorithm, named LSPR, has similar procedure with the ARSER algorithm. LSPR combines the spectral estimation and curve fitting techniques. For a given irregularly sampled time series, our method first estimates the period by Lomb-Scargle periodogram in frequency domain, and then models the periodic signals by harmonic regression method in time domain. Such joint strategy overcomes the limitations of Lomb-Scargle periodogram, and gives better descriptions for periodic patterns. Experiments on well-defined simulation data and \emph{Arabidopsis} and Yeast expression data show LSPR is more accurate in detecting periodicity than Lomb-Scargle periodogram and COSOPT for unevenly sampled short time-course expression data.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Results and Discussion %%
%%
\section*{Results and Discussion}
  \subsection*{Simulated data}
  \subsubsection*{Published data set}
    In this study, we analyzed the temporal microarray data to search genes that are periodically expressed. We developed the LSPR algorithm based on the Lomb-Scargle periodogram. Lomb-Scargle periodogram is a periodicity detection technology proposed by \cite{glynn2006detecting} to analyze gene expression data, especially for unevenly sampled time-series. The authors applied Lomb-Scargle periodogram to a given set of  simulated time-series and concluded an estimated relationship between the number of time-points, $N$, and the statistical significance of periodicity, $p$-value, as follows:
    \begin{equation}
    \label{eq:01}
    N \approx 5[1-\log_{10}p]
    \end{equation}
    According to Equation (\ref{eq:01}), for a statistical significant rhythmic gene with the $p$-value threshold 0.05, the number of time-points of the expression profiles will need to be more than 20.  The diurnal and circadian time course studies are usually designed to collect data every 4 hours over a course of 48 hours, generating expression profiles with 12 or 13 time-points\cite{Yamada2007}. Such short time-series make Lomb-Scargle periodogram fail to detect the periodic trends.

    We firstly compared our algorithm with COSOPT and Lomb-Scargle periodogram by testing them on a group of short periodic time-series (Table 1). This data set was generated by \cite{Michael2008} with five periodic patterns based on studies available in the literature. This data set contains 120 time series in all, with 24 samples for each periodic pattern respectively. We found Lomb-Scargle identified none of the periodic signals under the false discovery rate (FDR) threshold 0.05. The testing time-series contained 12 time-points, which led the statistical $p$-value of Lomb-Scargle periodogram not significant for such short time-series. COSOPT detected the periodicity by curve fitting, which processed more power than Lomb-Scargle method. It identified ~72\% samples from all the periodic signals.

    LSPR employed harmonic regression to estimate the $p$-value, which is a different way from Lomb-Scargle periodogram. Similar estimation between $N$ and $p$-value was conducted by using LSPR, as follows:
    \begin{equation}
    \label{eq:02}
    N \approx 5.6-\log_{10}p
    \end{equation}
    According to Equation (\ref{eq:02}), LSPR can identify a periodically expressed gene as significant at p-value 0.05 from a time-series with only 7 time-points. We found LSPR identified the most periodic signals in all the three algorithms, especially identifying the most periodic signals with spike and box-like periodic patterns.

  \subsubsection*{Our generated data sets}
    To test the LSPR algorithm and compare its performance with prior methods, we prepared comprehensive data sets containing periodic and non-periodic samples. The periodic samples were generated by two models. One is the stationary periodic model (see method section for details) defined by cosine curve with constant amplitude and mean value. The stationary model was widely used to generate synthetic data in prior studies \cite{ptitsyn2006permutation, liew2007spectral, wichert2004identifying}. Considering the dampening effect of circadian rhythm \cite{10.1371/journal.pcbi.1000580},  the periodic data were also generated by another non-stationary model defined by cosine curve with exponentially damped amplitude and mean value. Comparing the stationary model, the non-stationary one is more likely to approximate the natural biological rhythm \cite{refinetti2004non}.

    The non-periodic data in the synthetic data sets were generated by two random processes. One is white noise following standard normal distribution. Another widely accepted stochastic process in time series is autoregressive model considering the general correlation between successive measurements. Here we used autoregressive process of order one (AR(1)) to generate the random time-series.

    In our synthetic data sets, the periodic signals include 10 000 stationary and 10 000 non-stationary samples with varied period, phase and signal-to-noise ratio (SNR). The non-periodic signals include 10 000 white noise signals and 10 000 AR(1) signals. Each time-series possessed 20 unevenly sampled time-points over the course of two days. This sampling design is consist with the way taken by \cite{smith2004diurnal} (see method section for details).

    Four data sets were constructed by combining the periodic and non-periodic signals. We applied LSPR, Lomb-Scargle periodogram and COSOPT to analyze these four data sets. Since the task is to separate periodic signals from non-periodic signals, we can treat it as a binary classification problem. This allows receive operating characteristic (ROC) analysis to take place. Figure (1A-D) showed the sensitivity vs. (1 - specificity) for three algorithms as their determination thresholds (FDR $q$-value for LSPR and Lomb-Scargle, pMMC-$\beta$ for COSOPT) are varied. We found LSPR performed the best of the three algorithms in all the cases.

  \subsection*{Molecular data}
  \subsubsection*{\emph{Arabidopsis} diurnal expression data}
    Periodic processes (such as circadian or diurnal cycles) play an important role in regulating plant metabolic and physiological processes. Here two independent data sets were selected from the previous studies of diurnal gene expression of model plant, \emph{Arabidopsis}. These data sets were named after their respective first author: Smith data\cite{smith2004diurnal} and Blasing data\cite{blasing2005sugars}. Smith data were generated by harvesting RNA samples at eleven different time-points as follows: 0, 1, 2, 4, 8, 12, 13, 14, 16, 20, and 24 h (where time 0 is the onset of dark and 12 h is the onset of light). The 24 h time point is a repeat of 0 h. There are two biological replicates for each sample. Blasing data were generated by harvesting RNA samples at 6 time-points every 4 hours beginning at the end of the night with 3 biological replicates. These two data sets have similar experimental conditions. Application of Lomb-Scargle periodogram to the expression data identified zero transcripts as rhythmic in both Smith and Blasing data at the $q$-value cutoff 0.05. It appears that Lomb-Scargle periodogram has weak statistical power with such short signal length. Thus, we applied COSOPT and LSPR to identify rhythmic transcripts from Smith and Blasing data. Of all 22 810 \emph{Arabidopsis} genes, we found LSPR identified 7815 transcripts (35\% of \emph{Arabidopsis} whole genome) as rhythmic genes in the Blasing data and 6709 transcripts (30\% of \emph{Arabidopsis} whole genome) in the Smith data with $q$-value cutoff 0.05 (figure 2A). These fractions of clock-regulated genes were consistent with an estimation of between 31\% and 41\% of expressed genes being circadian-regulated reported by a recent study \cite{covington2008global}. The overlap of identified rhythmically expressed genes between Smith and Blasing data is 19\% of all expressed genes (figure 2A), which is higher than that of COSOPT's (figure 2B).

    Remarkably, the LSPR method covered 99\% and 97\% of rhythmic genes identified by COSOPT in Blasing and Smith data respectively. In addition, LSPR newly identified 9\% and 7\% of whole genome transcripts as rhythmic genes in Blasing and Smith data (figure 2C and 2D).

    \subsubsection*{Yeast cell cycle data}
    Generally, periodicity detection algorithms are applied to circadian rhythm and cell cycle expression data. We have tested the performance of our method in \emph{Arabidopsis} diurnal expression data. We will continue to validate its power in analyzing Yeast cell cycle expression data generated by Spellman \emph{et al}. \cite{spellman1998comprehensive}. In our benchmark, we included three computational methods: LSPR, Lomb-Scargle periodogram and the spectral estimation method proposed by Alan \emph{et al}. \cite{liew2007spectral}. As suggested by the study of Lichtenberg \emph{et al}. \cite{de2005comparison}, we provide the benchmark of these methods by applying them to three Yeast experimental data sets: alpha, cdc15 and cdc28. For Alpha and cdc15 experiments are raw data as obtained from \url{http://genome-www.stanford.edu/cellcycle/}. For cdc28 experiment are renormalized data generated by Lichtenberg \emph{et al}. \cite{de2005comparison}. Alpha experiment harvested samples at 7 min intervals for 119 min with a total of 18 time-points. Cdc15 data were measured every 10 min for 290 min, lacking observations for 20, 40, 60, 260 and 280 min time-points, and gives a total 24 time-points. For cdc28 experiment, samples were taken every 10 min for 160 min with a total of 17 time-points. Alpha and cdc28 were treated as evenly sampled data, while cdc15 were unevenly sampled data. In our analysis, we removed the genes with missing values for all sample points and obtained 4489 genes for alpha, 4381 genes for cdc15 and 6214 genes for cdc28.

    To measure the performance of different algorithms, Lichtenberg \emph{et al}. \cite{de2005comparison} proposed three benchmark sets B1, B2 and B3. B1 contains a total of 113 genes previously identified as periodically expressed in small-scale experiments. B2 contains 352 genes whose promoters were bound by at least one of the 9 known cell cycle transcription factors in two Chromatin IP studies, and therefore many of the genes in this benchmark set should be expected to be cell cycle regulated. B3 contains 518 genes annotated in MIPS \cite{Mewes01012002} as "cell cycle and DNA processing". However, since a large number of genes involved in the cell cycle are not subjected to transcriptional regulation (not periodic) and genes found in B1 were explicitly removed, only a small fraction of the genes in B3 are expected to be periodically expressed.

    Figure 3 shows the performance of each method on three Yeast cell cycle data sets. There is no single method outperforms than the others across all benchmark sets in all experiments. Apparently, all methods perform significantly better than random. LSPR and Lomb-Scargle methods give similar performance for cell cycle data, which is different from the results of analyzing circadian expression data. This demonstrates that Lomb-Scargle periodogram is more efficient to analyze cell cycle expression data other than circadian expression data, since the cell cycle experiments have more sampling time points. According to the comparisons of three methods on cdc15 and cdc28, both LSPR and Lomb-Scargle methods give better performance than the method proposed by Alan and his co-workers.

%%%%%%%%%%%%%%%%%%%%%%
\section*{Conclusions}
  We have designed the LSPR algorithm for periodicity detection, which is a three-step algorithm based on Lomb-Scargle periodogram and harmonic regression. LSPR features in identifying periodically expressed genes from unevenly sampled time-course expression profiles with short length. When comparing LSPR with well-established algorithms in well-defined synthetic data and microarray data for diurnal and cell cycle experiments, the results illustrated that LSPR was superior in identifying periodic patterns as well as robustness to short time-series.



%%%%%%%%%%%%%%%%%%
\section*{Methods}
  \subsection*{Overview of LSPR algorithm}
    LSPR is a three-step integrated algorithm, featuring in high resolution for short time-series and precise period detection for unevenly sampled data. Gene expression profiles are inevitably blended with formidable confounds or consist of linear trends. To remedy the inaccuracy causing by noise and linear trend, LSPR first employs two preprocessing techniques - trend removal and noise filtering prior to identifying period components. Subsequently, LSPR estimates the spectrum in frequency domain by Lomb-Scargle periodogram, which enable itself to extract periodic components from unevenly sampled time-series. Finally, for purpose of validating the statistical significance of identified periods, LSPR models the cyclic components of time-series by harmonic regression. Harmonic regression models involve estimating the amplitude, phase and mean value that best fit the time series. By following above three steps, periodic time-series could be discerned and LSPR uses false discovery rate method to deal with multiple testing problem for large-scale temporal microarray data.

  \subsection*{Data preprocessing strategies}
     Any linear trend in the time series has been removed by fitting the data to a linear regression curve, then subtracting this linear component from each time point.  Noise filtering is performed by smoothing the detrended time series under 4th degree polynomial procedure with the application of Savitzky-Golay filtering algorithm \cite{savitzky1964smoothing}, which is a low-pass filter that can efficiently remove pseudo-peaks in a spectrum caused by noise.

  \subsection*{Period detection}
  \subsubsection*{Lomb-Scargle Periodogram}
    Our method employs Lomb-Scargle periodogram to estimate the spectrum in frequency domain for unevenly spaced time series. This method was developed by Lomb and additionally elaborated by Scargle \cite{william1988numerical}, which makes use of each time point other than computing based on even intervals. For a time series $x_{i}\equiv x(t_{i})(i=1,2,...,N)$, where $t_{i}$ is the observation number. Lomb-Scargle periodogram implements an estimate of spectral power by
    \begin{equation}
    \label{eq:03}
    P_x(\omega) = \frac{1}{2}
\left(
  \frac { \left[ \sum_i x_i \cos \omega ( t_i - \tau ) \right] ^ 2}
        { \sum_i \cos^2 \omega ( t_i - \tau ) }
+
 \frac {\left[ \sum_i x_i \sin \omega ( t_i - \tau ) \right] ^ 2}
        { \sum_i \sin^2 \omega ( t_i - \tau ) }
\right)
    \end{equation}
    where $\bar{x}$ is mean; $\sigma^{2}$ is variance; $\omega$ is corresponding frequency in radians ($\omega=2\pi f$) and the constant $\tau$ is defined by the relation:
    \begin{equation}
    \label{eq:04}
    \tan(2\omega\tau)=\frac{\Sigma_{i}\sin 2\omega t_{i}}{\Sigma_{i}\cos 2\omega t_{i}}
    \end{equation}
    The generated periodogram from a cyclic gene expression profile will exhibit one or several significant peaks accounting for the periodic components. Any frequency $f$ corresponding to the peak is selected from a defined period window and transformed to time domain as estimated period (if any) by the formula $T=1/f$.

  \subsubsection*{Procedure of period detection}
    LSPR obtains the period by using the following step-by-step procedure:
\begin{enumerate}
 \item Remove the linear trend in time-series $\{x_{i}\}$, denoting the detrended time-series as $\{\dot{x}_{i}\}$.
 \item Smooth $\{\dot{x}_{i}\}$ by a fourth-order Savitzky-Golay algorithm. The smoothed time-series is denoted as $\{\ddot{x}_{i}\}$.
 \item Calculate the Lomb-Scargle periodogram of $\{\ddot{x}_{i}\}$ by Equation~(\ref{eq:03}), and select all periods $\{\ddot{T}_{j}\}\in[20,28]$ that show peaks in the spectrum.
 \item Calculate the Lomb-Scargle periodogram of detrended time-series $\{\dot{x}_{i}\}$ by Equation~(\ref{eq:03}), and select all periods $\{\dot{T}_{j}\}\in[20,28]$ that show peaks in the spectrum.
 \item The periods $\{\dot{T}_{j}\}$ and $\{\ddot{T}_{j}\}$ are chosen as input to the harmonic regression analysis for $\{\dot{x}_{i}\}$ by Equation~(\ref{eq:05}).

 \item The optimum periods in $\{x_{i}\}$ are determined by Akaike's information criterion \cite{akaike1974new} among the regression models generated in step 5.
 \end{enumerate}

\subsection*{Harmonic regression}
  After periodic components being identified, harmonic regression models are used to represent the cyclic trends by fitting the time-series with sinusoidal functions as following:
  \begin{equation}
  \label{eq:05}
    x_{i} = \mu + \sum_{j=1}^{K}\beta_{j}\cos(\frac{2\pi}{T_{j}}t_{i} + \varphi_{j}) + \varepsilon_{i}
  \end{equation}
    Where $x_{i}$ is observed value time $t_{i}$; $\mu$ is the mean level of the time series, $\beta_{j}$ are amplitude of different sinusoidal waveforms;  $\varphi_{j}$ are the phase or the location of peaks relative to time zero; $\varepsilon_{i}$ is uncorrelated random variable; $t_{i}$ are sampling time-points, and $T_{j}$ are estimated periods.

    Equation~\ref{eq:05} can be reduced to a simple linear equation:
    \begin{equation}
    \label{eq:06}
    x_{i} = \mu + \sum_{j=1}^{K}\{p_{j}\cos(\frac{2\pi}{T_{j}}t_{i}) + q_{j}\sin(\frac{2\pi}{T_{j}}t_{i})\} + \varepsilon_{i}
    \end{equation}
    where $p_{j} = \beta_{j}\cos\varphi_{j}$, $q_{j} = -\beta_{j}\sin\varphi_{j}$, ordinary least squares (OLS) regression procedure is used to estimate parameters $p_{j}$, $q_{j}$, $\mu$; Amplitude $\beta_{j} = \sqrt{p_{j}^{2} + q_{j}^2}$ and phase $\varphi$ can be obtained by $\tan\varphi_{j} = -q_{j}/p_{j}$.

    By applying the harmonic regression, periodicity can be fully described by four parameters: period, phase, amplitude and mean level. F-test for $\beta_{j}$ is employed to computing the statistical significance of the model.

 \subsection*{Multiple testing correction}
  To analyze the large-scale temporal expression profiles, LSPR uses false discovery rate method proposed by \cite{Storey2003} for multiple testing correction. Briefly, by examining the distribution of p-values from the given data set, an estimate of the proportion that is truly non-rhythmic can be derived. The p-value for each transcript can be converted to a more stringent q-value which represents the false discovery rate. In our study, we consider genes with q-value $<$ 0.05 to be rhythmically expressed.

 \subsection*{Simulation data}
 In our numerical experiment part, we prepared periodic and non-periodic data sets to test the effectiveness of LSPR. To generate the periodic time-series, we applied two models to simulate the expression data regulated by periodic genes. The stationary model to represent an ``ideal'' gene expression is defined as:
 \begin{equation}
\label{eq:07}
x_{i} = \mathrm{SNR}\cdot2\cos(\frac{2\pi}{T}t_{i} - \varphi) + \varepsilon_{i}
\end{equation}
where $\mathrm{SNR}$ is signal-to-noise ratio; $\mathrm{T}$ is period; $\varphi$ is phase; and $\varepsilon_{i}$ is $(\mu=0,\sigma=1)$ normally distributed noise terms.

The other is non-stationary model which considers the dampening effect of free-running rhythms. It is defined as:
\begin{equation}
\label{eq:08}
x_{i} = 500\cdot e^{-0.01t_{i}} + \mathrm{SNR}\cdot100\cdot e^{-0.01t_{i}}\cdot\cos(\frac{2\pi}{T}t_{i} - \varphi) + \varepsilon_{i}
\end{equation}
where $\varepsilon_{i}$ is $(\mu=0,\sigma=50)$ normally distributed noise; the mean level and amplitude exponentially decay over time. $t_{i}$ are unevenly sampled at 0h, 1h, 2h, 4h, 8h, 12h, 13h, 14h, 16h, 20h with one replicate for each time point respectively. By varying $\mathrm{SNR}$, $T$ and $\varphi$ in above two models, multiple periodic time-series are generated.
Non-periodic signals are generated from two random processes: standard normal distributed white noise and AR(1) model defined by:
\begin{equation}
\label{eq:09}
x_{t} = c + \alpha x_{t-1}+ \varepsilon_{t}
\end{equation}
where $\alpha$ is AR coefficient, $c$ is a constant and $\varepsilon_{t}$ is white noise.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Authors contributions}
    Rendong Yang conceived the experiments and prepared the data for testing and analyzing. Chen Zhang developed the LSPR program and performed data analysis. Rendong Yang and Chen Zhang tested the algorithm and wrote the paper. Zhen Su supervised the whole project from algorithm development, data preparation to drafting the manuscript.



%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Acknowledgements}
  \ifthenelse{\boolean{publ}}{\small}{}
  We thank Alan Wee-Chung Liew for sharing codes in his study; Daofeng Li for web server assistance,  Wenying Xu, Yi Ling for helpful advices and discussions.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                  The Bibliography                       %%
%%                                                         %%
%%  Bmc_article.bst  will be used to                       %%
%%  create a .BBL file for submission, which includes      %%
%%  XML structured for BMC.                                %%
%%  After submission of the .TEX file,                     %%
%%  you will be prompted to submit your .BBL file.         %%
%%                                                         %%
%%                                                         %%
%%  Note that the displayed Bibliography will not          %%
%%  necessarily be rendered by Latex exactly as specified  %%
%%  in the online Instructions for Authors.                %%
%%                                                         %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


{\ifthenelse{\boolean{publ}}{\footnotesize}{\small}
 \bibliographystyle{bmc_article}  % Style BST file
  \bibliography{bmc_article} }     % Bibliography file (usually '*.bib' )

%%%%%%%%%%%

\ifthenelse{\boolean{publ}}{\end{multicols}}{}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Figures                       %%
%%                               %%
%% NB: this is for captions and  %%
%% Titles. All graphics must be  %%
%% submitted separately and NOT  %%
%% included in the Tex document  %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
%% Do not use \listoffigures as most will included as separate files

\section*{Figures}

  \subsection*{Figure 1 - ROC analysis of LSPR, COSOPT and Lomb-Scargle periodogram}
  Testing data sets contains: (A) 10000 stationary periodic signals and 10000 white noise, (B) 10000 non-stationary periodic signals and 10000 white noise, (C) 10000 stationary periodic signals and 10000 AR(1)-based signals and (D) 10000 non-stationary periodic signals and 10000 AR(1)-based signals. Greater area under ROC curve means more accuracy for the prediction of the algorithm. Relatively, LSPR keeps the lowest false positive rate and false negative rate in all the cases.

  \subsection*{Figure 2 - Comparisons of periodicity detection algorithms in two \emph{Arabidopsis} data sets.}
  (A) Overlap of periodic transcripts identified by LSPR for Blasing and Smith data. (B) Overlap of periodic transcripts identified by COSOPT for Blasing and Smith data. (C) Overlap of periodic transcripts in Blasing data identified by LSPR and COSOPT. (D) Overlap of periodic transcripts in Smith data identified by LSPR and COSOPT. Comparing with COSOPT, LSPR is more sensitive to unevenly spaced data (16\% vs. 11\%) and LSPR detected more periodic genes (9\% in Blasing data and 7\% in Smith data).

  \subsection*{Figure 3 - Comparison of periodicity detection algorithms in Yeast cell cycle data}
  The fraction of benchmark genes contained in the top 800 ranked list is plotted for each algorithm, benchmark set (B1, B2, B3) and experiment (alpha, cdc15, cdc28), respectively. The methods are colored as follows: LSPR (our method, orange), Lomb-Scargle (\cite{glynn2006detecting}, magenta) and Alan \emph{et al}. (\cite{liew2007spectral}, green). Each ranked list is formed by sorting the p-value of gene profiles in the ascending order. Random performance could be observed in the black dot.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Tables                        %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Use of \listoftables is discouraged.
%%
\section*{Tables}
  \subsection*{Table 1 - Number of inferred periodic samples from short time-series with multiple periodic patterns}
    \par
    \mbox{
      \begin{tabular}{|c|c|c|c|c|c|c|}
        \hline Method & & & Waveforms & & & \\
        \hline  & spike (24) & box1 (24) & box2 (24) & rigid (24) & sine (24) & random (24) \\ \hline
        Lomb-Scargle & 0  & 0 & 0  & 0 & 0 & 0 \\ \hline
        COSOPT & 7 & 24 & 7 & 24 & 24 & 1 \\ \hline
        LSPR & 16  & 24 & 16  & 24 & 24 & 0 \\ \hline
      \end{tabular}
      }

\end{bmcformat}
\end{document}







